/*
 * Copyright (c) 2011 Seoul National University
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are
 * met: redistributions of source code must retain the above copyright
 * notice, this list of conditions and the following disclaimer;
 * redistributions in binary form must reproduce the above copyright
 * notice, this list of conditions and the following disclaimer in the
 * documentation and/or other materials provided with the distribution;
 * neither the name of the copyright holders nor the names of its
 * contributors may be used to endorse or promote products derived from
 * this software without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
 * OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 *
 * Authors: Jung Ho Ahn  ( gajh@snu.ac.kr )
 *
 */

#include <cassert>
#include <cstring>
#include <pthread.h>
#include <sstream>
#include <vector>
#include "SequenceBuffer.h"
#include "KmerHashTable.h"
#include "Kmer.h"
#include "Reader.h"
#include "FastAReader.h"
#include "FastQReader.h"
#include "BarrierUtil.h"

using namespace std;

pthread_mutex_t read_file_m;
pthread_mutex_t hashmaps_m;
pthread_barrier_t global_barrier;

void * linkKmersParallel(void *);
void * populateKmerNodesParallel(void *);
void * processReadsParallel(void *);
void * initializeContigsParallel(void *);
void * initializeContigsParallelTight(void *);


ostream & operator<<(ostream & out, const KmerCoverageTable & k)
{
  uint64_t size = 0;

  out << "KmerCoverageTable: hashsize= " << k.getTableSize() 
      << ", tmphashsize= " << k.getTmpTableSize() << ", kmers=" << k.getKmerNumber();

  for (int i = 0; i < k.numThreads; i++)
  {
    size += k.standaloneContigs[i].size();
  }

  out << ", standaloneContigSize= "  << size;

  size = 0;
  for (int i = 0; i < k.numThreads; i++)
  {
    size += k.activeContigs[i].size();
  }

  out << ", activeContigSize= "      << size;
  return out;
}


KmerCoverageTable::KmerCoverageTable(InputOptions * options_)
 :options(options_), kcount(0), ccount(0), histogram()
{
  string contigFilename      = options->get_param_str("output_contig_filename");
  string outputFastaFilename = options->get_param_str("output_fasta_filename");
  ctgout.open(contigFilename.c_str());
  fastaout.open(outputFastaFilename.c_str());

  if (ctgout.good() == false)
  {
    cout << "cannot open " << contigFilename << " for storing contigs." << endl;
    exit(1);
  }

  if (fastaout.good() == false)
  {
    cout << "cannot open " << contigFilename << " for displaying output contigs." << endl;
    exit(1);
  }

  kmersize       = options->get_param_uint64("kmer_size", 29);
  maxErrorCount  = options->get_param_uint64("max_count_treated_as_sequencer_error", 1);
  minRepeatCount = options->get_param_uint64("min_count_treated_as_repeat", 10000);
  minValidContigLength = options->get_param_uint64("min_valid_contig_length", 100);
  createContigSizeHistogram = options->get_param_bool("create_contig_size_histogram");

  Kmer::lenQ     = kmersize / nNtPerUnit;
  Kmer::lenR     = kmersize % nNtPerUnit;

  uint32_t bufSize   = options->get_param_uint64("seq_buffer_size",     25000);
  uint32_t maxSeqLen = options->get_param_uint64("max_sequence_length", 256);

  display_repeats             = options->get_param_bool("display_repeats");
  per_file_error_removal      = options->get_param_bool("per_file_error_removal");
  per_size_error_removal      = options->get_param_bool("per_size_error_removal");
  error_removal_block_size    = options->get_param_uint64("error_removal_block_size", 100000);
  error_removal_block_size_first = options->get_param_uint64("error_removal_block_size_first", error_removal_block_size);
  first_error_removal         = true;

  parallelize_kmer_linking    = options->get_param_bool("parallelize_kmer_linking", true);
  parallelize_kmer_node_population = options->get_param_bool("parallelize_kmer_node_population", true);
  parallelize_read_processing = options->get_param_bool("parallelize_read_processing", true);
  parallelize_contig_initialization = options->get_param_bool("parallelize_contig_initialization", true);
  tight_memory_management     = options->get_param_bool("tight_memory_management");
  display_fasta               = options->get_param_bool("display_fasta", true);
  display_confusing_kmers     = options->get_param_bool("display_confusing_kmers");
  disconnect_confusing_kmers  = options->get_param_bool("disconnect_confusing_kmers");
  kmerBatchSize               = options->get_param_uint64("kmer_batch_size", 1048576);
  interactive_mode            = options->get_param_bool("interactive_mode");
  minimal_hitcount_for_confusing_kmer = options->get_param_uint64("minimal_hitcount_for_confusing_kmer", 20);

  numHashmaps        = options->get_param_uint64("num_hashmaps",        16384);
  numThreads         = options->get_param_uint64("num_threads",             1);
  dotPerLine         = options->get_param_uint64("dot_per_line",          100);
  dotScaling         = options->get_param_uint64("dot_scaling",           100);
  istringstream dNThreshold(options->get_param_str("dominant_node_threshold"));
  BlockAllocatorMT::block_size = options->get_param_uint64("block_allocator_size", 33554432);  // 32MB

  if (dNThreshold.good() == true)
  {
    dNThreshold >> dominantNodeThreshold;
  }
  else
  {
    dominantNodeThreshold = 0.25;
  }

  //bamt.initiate();
  hashmapstmp        = NULL;
  hashmaps           = new KmerHash[numHashmaps];
  hashmap_m          = new pthread_mutex_t[numHashmaps];
  graphs             = new list<KmerNode *>[numHashmaps];

  standaloneContigs  = new ContigHash[numThreads];
  activeContigs      = new ContigHash[numThreads];

  for (int i = 0; i < numHashmaps; i++)
  {
    pthread_mutex_init(&(hashmap_m[i]), NULL);
  }

  pthread_mutex_init(&read_file_m, NULL);
  pthread_mutex_init(&hashmaps_m, NULL);
  pthread_barrier_init(&global_barrier, NULL, numThreads);

  seqbuf = new SequenceBuffer *[numThreads];
  blockAllocators = new BlockAllocatorST *[numThreads];
  for (int i = 0; i < numThreads; i++)
  {
    seqbuf[i] = new SequenceBuffer(bufSize, maxSeqLen);
    blockAllocators[i] = new BlockAllocatorST(BlockAllocatorMT::block_size);
  }

  if (parallelize_kmer_node_population == true ||
      parallelize_kmer_linking == true ||
      parallelize_read_processing == true ||
      parallelize_contig_initialization)
  {
    threads = new pthread_t[numThreads];
    th_args = new kmer_thread_args[numThreads];

    for (int i = 0; i < numThreads; i++)
    {
      th_args[i].table = this;
      th_args[i].tid   = i;
    }
  }
}


KmerCoverageTable::~KmerCoverageTable()
{
  for (int i = 0; i < numHashmaps; i++)
  {
    pthread_mutex_destroy(&(hashmap_m[i]));
  }

  pthread_mutex_destroy(&read_file_m);
  pthread_mutex_destroy(&hashmaps_m);

  if (parallelize_kmer_node_population == true ||
      parallelize_kmer_linking == true ||
      parallelize_read_processing == true)
  {
    delete [] threads;
    delete [] th_args;
  }
}


void KmerCoverageTable::processARead(const char * sequence)
{
  uint32_t len = strlen(sequence);
  int      lastInvalid = kmersize;
  //char     latestChar  = sequence[0];
  //uint32_t numSameChar = 0;

  if (len < kmersize)
  {
    return;	// no kmer can be generated
  }

  Kmer kmer;
  Kmer kmerRC;  // RC stands for reverse complement

  // Initialize kmer
  for (uint32_t i = 0; i < len; i++)
  {
    lastInvalid = (Kmer::isValidNucleotide(sequence[i]) == false) ? kmersize : lastInvalid - 1;
    kmer.pushFromLSB(Kmer::charToNucleotide(sequence[i]));
    kmerRC.pushFromMSB(Kmer::charToComplementaryNucleotide(sequence[i]));

    //numSameChar = (latestChar == sequence[i]) ? (numSameChar + 1) : 1;
    //latestChar  = sequence[i];

    if (lastInvalid <= 0)// && numSameChar < kmersize)
    {
      add(&kmer, &kmerRC);
    }
  }
}


void KmerCoverageTable::add(Kmer * kmer, Kmer * kmerRC)
{
  if (*kmerRC < *kmer)
  {
    Kmer * temp = kmerRC;
    kmerRC = kmer;
    kmer   = temp;
  }

  int hashValue         = kmer->getLSW() % numHashmaps;
  KmerHash & hashmap    = hashmaps[hashValue];
  KmerHash::iterator it = hashmap.find(kmer);

  if (it == hashmap.end())
  {
    //Kmer * newkmer     = new Kmer(*kmer);
    Kmer * newkmer = (Kmer *)blockAllocators[0]->allocate(sizeof(Kmer));
#ifdef __GENERAL_KMER__
    newkmer->kmer = (KmerUnit *)blockAllocators[0]->allocate(sizeof(KmerUnit)*(Kmer::lenQ+1));
#endif
    *newkmer = *kmer;
    hashmap.insert(KmerPair(newkmer, CountOrNodes()));
  }
  else
  {
    ((it->second).count)++;
  }
}


void KmerCoverageTable::addParallel(Kmer * kmer, Kmer * kmerRC, int tid)
{
  if (*kmerRC < *kmer)
  {
    Kmer * temp = kmerRC;
    kmerRC = kmer;
    kmer   = temp;
  }

  int hashValue = kmer->getLSW() % numHashmaps;
  KmerHash::iterator it;

  pthread_mutex_lock(&(hashmap_m[hashValue]));
  KmerHash & hashmap = hashmaps[hashValue];
  it = hashmap.find(kmer);

  if (it == hashmap.end())
  {
    //Kmer * newkmer     = new Kmer(*kmer);
    Kmer * newkmer = (Kmer *)blockAllocators[tid]->allocate(sizeof(Kmer));
#ifdef __GENERAL_KMER__
    newkmer->kmer = (KmerUnit *)blockAllocators[tid]->allocate(sizeof(KmerUnit)*(Kmer::lenQ+1));
#endif
    *newkmer = *kmer;
    hashmap.insert(KmerPair(newkmer, CountOrNodes()));
  }
  else
  {
    ((it->second).count)++;
  }
  pthread_mutex_unlock(&(hashmap_m[hashValue]));
}


void KmerCoverageTable::addParallelPrune(Kmer * kmer, Kmer * kmerRC, int tid)
{
  if (*kmerRC < *kmer)
  {
    Kmer * temp = kmerRC;
    kmerRC = kmer;
    kmer   = temp;
  }

  int hashValue = kmer->getLSW() % numHashmaps;
  KmerHash::iterator it;

  pthread_mutex_lock(&(hashmap_m[hashValue]));
  KmerHash & hashmap = hashmaps[hashValue];
  it = hashmap.find(kmer);

  if (it == hashmap.end())
  {
    KmerHashTmp * hashmaptmp = hashmapstmp[hashValue];
    KmerHashTmp::iterator ittmp = hashmaptmp->find(*kmer);

    if (ittmp == hashmaptmp->end())
    {
      hashmaptmp->insert(KmerPairTmp(*kmer, 1));
    }
    else if (ittmp->second == maxErrorCount)
    {
      //Kmer * newkmer     = new Kmer(*kmer);
      Kmer * newkmer = (Kmer *)blockAllocators[tid]->allocate(sizeof(Kmer));
#ifdef __GENERAL_KMER__
      newkmer->kmer = (KmerUnit *)blockAllocators[tid]->allocate(sizeof(KmerUnit)*(Kmer::lenQ+1));
#endif
      *newkmer = *kmer;
      hashmap.insert(KmerPair(newkmer, CountOrNodes(maxErrorCount+1)));
    }
    else
    {
      (ittmp->second)++;
    }
  }
  else
  {
    ((it->second).count)++;
  }
  pthread_mutex_unlock(&(hashmap_m[hashValue]));
}


KmerNode * KmerCoverageTable::find(Kmer * kmer) const
{
  KmerHash::const_iterator it;
  KmerHash & hashmap = hashmaps[kmer->getLSW() % numHashmaps];
  it = hashmap.find(kmer);
  if (it == hashmap.end())
  {
	return NULL;
  }
  else
  {
    return it->second.nodes;
  }
}


void KmerCoverageTable::initializeContigs()
{
  if (parallelize_contig_initialization == true)
  {
    nDots         = 0;
    countHashmaps = 0;
    for (int i = 0; i < numThreads; i++)
    {
      if (tight_memory_management == true)
      {
        pthread_create(&(threads[i]), NULL, initializeContigsParallelTight, (void *)(&(th_args[i])));
      }
      else
      {
        pthread_create(&(threads[i]), NULL, initializeContigsParallel, (void *)(&(th_args[i])));
      }
    }

    for (int i = 0; i < numThreads; i++)
    {
      pthread_join(threads[i], NULL);
    }

    cout << endl;
    ctgout << endl;
    return;
  }

  for (int i = 0; i < numHashmaps; i++)
  {
    KmerHash & hashmap = hashmaps[i];
    KmerHash::iterator iter = hashmap.begin();

    while (iter != hashmap.end())
    {
      KmerNode * pnode = (iter->second).nodes;

      // a node with a branch is a graph node
      // if pnode->kmer == NULL, kmer is moved to the ContigNode already. simply delete it
      if (pnode->hasBranch() == false && pnode->valid == true)
      {
        uint64_t pointer = blockAllocators[0]->allocate(sizeof(ContigNode));
        ContigNode * contigNode = new(pointer) ContigNode(pnode, ccount++);
        printAndUpdateContigSizeHistogram(contigNode, 0);
      }

      pnode = (iter->second).nodes + 1;

      if (pnode->hasBranch() == false && pnode->valid == true)
      {
        uint64_t pointer = blockAllocators[0]->allocate(sizeof(ContigNode));
        ContigNode * contigNode = new(pointer) ContigNode(pnode, ccount++);
        printAndUpdateContigSizeHistogram(contigNode, 0, true);
      }
      ++iter;
    }
  }

  ctgout << endl;
}


void KmerCoverageTable::displayNodes()
{
  ctgout << "standalone contigs" << endl;
  for (int i = 0; i < numThreads; i++)
  {
    ContigHash & chash = standaloneContigs[i];
    for (ContigHash::iterator iter = chash.begin(); iter != chash.end(); ++iter)
    {
      ctgout << *(*iter) << endl;
    }
  }
  ctgout << endl;

  ctgout << "active contigs" << endl;
  for (int i = 0; i < numThreads; i++)
  {
    ContigHash & chash = activeContigs[i];
    for (ContigHash::iterator iter = chash.begin(); iter != chash.end(); ++iter)
    {
      ctgout << *(*iter) << endl;
    }
  }
  ctgout << endl;

  ctgout << "kmer nodes" << endl;
  for (int i = 0; i < numHashmaps; i++)
  {
    list<KmerNode *>::iterator iter;
    list<KmerNode *> & graph = graphs[i];

    for (iter = graph.begin(); iter != graph.end(); ++iter)
    {
      ctgout << *(*iter) << endl;
    }
  }
  ctgout << endl;
}


void KmerCoverageTable::displayKmerCountHistogram()
{
  KmerHash::iterator it;

  for (int i = 0; i < numHashmaps; i++)
  {
    KmerHash & hashmap = hashmaps[i];
    for (it = hashmap.begin(); it != hashmap.end(); it++)
    {
      uint64_t count = (it->second).count;

      if (count)
      {
        if (histogram.find(count) == histogram.end())
        {
          histogram[count] = 1;
        }
        else
        {
          histogram[count]++;
        }
      }
    }
  }

  cout << "*** Kmer Count Histogram ***" << endl;
  long counts = 0;

  for (map<uint64_t, uint64_t>::iterator iter = histogram.begin(); iter != histogram.end(); ++iter)
  {
    cout.width(5);
    cout << right << iter->first;
    cout << " : " << iter->second << endl;
    counts += iter->second;
  }
  cout << "Total kmer count: " << counts << endl;
}


void KmerCoverageTable::populateKmerNodes()
{
  if (display_repeats == true)
  {
    ctgout << "following kmers are ignored since their hitcounts >= " << minRepeatCount
           << " or <= " << maxErrorCount <<  endl;
  }

  if (parallelize_kmer_node_population == true)
  {
    nDots         = 0;
    countHashmaps = 0;
    for (int i = 0; i < numThreads; i++)
    {
      pthread_create(&(threads[i]), NULL, populateKmerNodesParallel, (void *)(&(th_args[i])));
    }

    for (int i = 0; i < numThreads; i++)
    {
      pthread_join(threads[i], NULL);
    }

    cout << endl;
    return;
  }

  uint64_t mapIdx = 0;
  const uint64_t mapPerDot  = 100000;

  // trim repeats and noises first
  if (display_repeats == true)
  {
    ctgout << "following kmers are ignored since their hitcounts >= " << minRepeatCount
           << " or <= " << maxErrorCount <<  endl;
  }

  for (int i = 0; i < numHashmaps; i++)
  {
    KmerHash & hashmap = hashmaps[i];
    KmerHash::iterator it = hashmap.begin();

    while (it != hashmap.end())
    {
      mapIdx++;
      if (mapIdx % mapPerDot == 0)
      {
        cout << "." << flush;

        if (mapIdx % (mapPerDot * dotPerLine) == 0)
        {
          cout << " populated " << mapIdx << " kmer nodes" << endl;
        }
      }

      assert(it->first != NULL && (it->second).count > 0);
      if ((it->second).count >= minRepeatCount || (it->second).count <= maxErrorCount)
      {
        if (display_repeats == true && (it->second).count >= minRepeatCount)
        {
          ctgout << *(it->first) << ", hitcount= " << (it->second).count << endl;
        }
        //delete it->first;
        hashmap.erase(it++);
      }
      else
      {
        uint64_t kmerNodePair = blockAllocators[0]->allocate(sizeof(KmerNode) * 2);
        uint64_t count = (it->second).count;
        /*KmerNode * newnode   = new KmerNode(it->first, kcount++, (it->second).count, false);
        KmerNode * newnodeRC = new KmerNode(it->first, kcount++, (it->second).count, true);
        // TODO: kcount could be a unnecessary serialization point
        it->second.nodes.kmerNode   = newnode;
        it->second.nodes.kmerNodeRC = newnodeRC;
        */
        it->second.nodes   = new(kmerNodePair)  KmerNode(it->first, kcount++, count, false);
        new(kmerNodePair+sizeof(KmerNode)) KmerNode(it->first, kcount++, count, true);
        ++it;
      }
    }
  }
  cout << " populated " << mapIdx << " kmers" << endl;
}


void KmerCoverageTable::linkKmers()
{
  countLinks = 0;
  nDots      = 0;

  if (parallelize_kmer_linking == true)
  {
    countHashmaps = 0;

    for (int i = 0; i < numThreads; i++)
    {
      pthread_create(&(threads[i]), NULL, linkKmersParallel, (void *)(&(th_args[i])));
    }

    for (int i = 0; i < numThreads; i++)
    {
      pthread_join(threads[i], NULL);
    }

    cout << endl << "OK, link count = " << countLinks << endl;
    return;
  }

  uint64_t mapIdx = 0;
  const uint64_t mapPerDot  = 100000;
  Kmer left, right, leftRC, rightRC;

  // link kmers
  for (int i = 0; i < numHashmaps; i++)
  {
    KmerHash::iterator it;
    KmerHash & hashmap = hashmaps[i];
    list<KmerNode *> & graph = graphs[i];

    for (it = hashmap.begin(); it != hashmap.end(); ++it)
    {
      mapIdx++;
      if (mapIdx % mapPerDot == 0)
      {
        cout << "." << flush;

        if (mapIdx % (mapPerDot * dotPerLine) == 0)
        {
          cout << " processed " << mapIdx << " kmers" << endl;
        }
      }

      KmerNode * kmerNode   = it->second.nodes;
      KmerNode * kmerRCNode = it->second.nodes + 1;
      assert(kmerNode != NULL && kmerRCNode != NULL);
      Kmer *     kmer       = it->first;
      //Kmer *     kmerRC     = kmerRCNode->getKmer();
      bool foundHasBranch = kmerNode->hasBranch();

      for (int j = 0; j < 4; j++)
      {
        right = *kmer;
        right.pushFromLSB((nucleotide_type)j);
        rightRC = *kmer;
        rightRC.reverseComplement();
        rightRC.pushFromMSB(Kmer::ComplementaryNucleotide((nucleotide_type)j));

        bool isRightRCSmall = (rightRC < right);
        KmerNode * rightNC = (isRightRCSmall) ? find(&rightRC) : find(&right);
        if (rightNC)
        {
          KmerNode * rightNode   = rightNC;
          KmerNode * rightRCNode = rightNC + 1;
          countLinks++;
          if (right == *(rightNode->kmer))
          {
            kmerNode->setNext(rightNode, (nucleotide_type)j);
            kmerRCNode->setPrev(rightRCNode, Kmer::ComplementaryNucleotide((nucleotide_type)j));
          }
          else
          {
            kmerNode->setNext(rightRCNode, (nucleotide_type)j);
            kmerRCNode->setPrev(rightNode, Kmer::ComplementaryNucleotide((nucleotide_type)j));
          }
          if (foundHasBranch == false && kmerNode->hasBranch() == true)
          {
            foundHasBranch = true;
            graph.push_back(kmerNode);
            graph.push_back(kmerRCNode);
          }
        }

        left  = *kmer;
        left.pushFromMSB((nucleotide_type)j);
        leftRC = *kmer;
        leftRC.reverseComplement();
        leftRC.pushFromLSB(Kmer::ComplementaryNucleotide((nucleotide_type)j));

        bool isLeftRCSmall = (leftRC < left);
        KmerNode * leftNC = (isLeftRCSmall) ? find(&leftRC) : find(&left);
        if (leftNC)
        {
          KmerNode * leftNode   = leftNC;
          KmerNode * leftRCNode = leftNC + 1;
          countLinks++;
          if (left == *(leftNode->kmer))
          {
            kmerNode->setPrev(leftNode, (nucleotide_type)j);
            kmerRCNode->setNext(leftRCNode, Kmer::ComplementaryNucleotide((nucleotide_type)j));
          }
          else
          {
            kmerNode->setPrev(leftRCNode, (nucleotide_type)j);
            kmerRCNode->setNext(leftNode, Kmer::ComplementaryNucleotide((nucleotide_type)j));
          }
          if (foundHasBranch == false && kmerNode->hasBranch() == true)
          {
            foundHasBranch = true;
            graph.push_back(kmerNode);
            graph.push_back(kmerRCNode);
          }
        }
      }
    }
  }
  cout << " processed " << mapIdx << " kmers" << endl;
  cout << "OK, link count = " << countLinks << endl;
}


void KmerCoverageTable::processReads()
{
  Reader * reader;
  string filename;

  if (per_file_error_removal == true && display_repeats == true)
  {
    ctgout << "following kmers are ignored since their hitcounts >= " << minRepeatCount
           << " or <= " << maxErrorCount <<  endl;
  }

  while (options->fasta_filenames.empty() == false || options->fastq_filenames.empty() == false)
  {
    if (options->fasta_filenames.empty() == false)
    {
      reader = new FastAReader();
      filename = options->fasta_filenames.front();
      options->fasta_filenames.pop_front();
    }
    else
    {
      reader = new FastQReader();
      filename = options->fastq_filenames.front();
      options->fastq_filenames.pop_front();
    }

    if (reader->openFile(filename.c_str()) == false)
    {
      cout << "failed to open file " << filename << endl;
      exit(1);
    }
    cout << endl << "Reading and processing " << filename << endl;

    if (parallelize_read_processing == true && 
        (per_file_error_removal == true || per_size_error_removal == true))
    {
      bamt.clearBlocks(true);
      bamt.initiate();
      hashmapstmp  = new KmerHashTmp * [numHashmaps];
      for (int i = 0; i < numHashmaps; i++)
      {
        uint64_t pointer = bamt.allocate(sizeof(KmerHashTmp));
        hashmapstmp[i] = new((void *)pointer) KmerHashTmp();
      }
    }


    countReads     = 0;
    countReadsCurr = 0;
    nDots          = 0;

    if (parallelize_read_processing == true)
    {
      for (int i = 0; i < numThreads; i++)
      {
        th_args[i].reader = reader;
        pthread_create(&(threads[i]), NULL, processReadsParallel, (void *)(&(th_args[i])));
      }

      for (int i = 0; i < numThreads; i++)
      {
        pthread_join(threads[i], NULL);
      }
    }
    else
    {
      int current = 0;

      while (true)
      {
        current = reader->readToBuffer(seqbuf[0]);

        if (current <= 0)
        {
          break;
        }
        else
        {
          countReads += current;
          cout << "." << flush;
          nDots++;
          if (nDots % dotPerLine == 0)
          {
            cout << " read " << countReads << " reads" << endl << flush;
          }
        }
        int lineCount = seqbuf[0]->getLineCount();
        for (int i = 0; i < lineCount; i++)
        {
          processARead(seqbuf[0]->getLinePointer(i));
        }
        seqbuf[0]->clear();
      }
    }

    cout << " read " << countReads << " reads" << endl << flush;
    cout << *this << endl;
    delete reader;
  }
  bamt.clearBlocks(false);
  bamt.initiate();
  hashmapstmp  = NULL;
}


uint64_t KmerCoverageTable::getTableSize() const
{
  uint64_t size = 0;

  for (int i = 0; i < numHashmaps; i++)
  {
    KmerHash & hashmap = hashmaps[i];
    size += hashmap.size();
  }

  return size;
}

uint64_t KmerCoverageTable::getKmerNumber() const
{
  uint64_t kmers = 0;
  for (int i=0; i<numHashmaps; i++)
  {
    KmerHash & hashMap = hashmaps[i];
    for (KmerHash::iterator it = hashMap.begin(); it != hashMap.end(); ++it) {
      kmers += (it->second).count;
    }
  }
  return kmers;
}


uint64_t KmerCoverageTable::getTmpTableSize() const
{
  if (hashmapstmp == NULL)
  {
    return 0;
  }

  uint64_t size = 0;

  for (int i = 0; i < numHashmaps; i++)
  {
    KmerHashTmp * hashmap = hashmapstmp[i];
    size += hashmap->size();
  }

  return size;
}


void KmerCoverageTable::displayContigSizeHistogram()
{
  if (createContigSizeHistogram == true)
  {
    cout << "*** Contig Size Histogram ***" << endl;
    uint64_t aggregateContigSize = 0;

    map<uint64_t, uint64_t>::reverse_iterator iter;
    for (iter = contigSizeHistogram.rbegin(); iter != contigSizeHistogram.rend(); ++iter)
    {
      cout << " + ";
      cout.width(5);
      cout << right << iter->first;
      cout << " * ";
      cout.width(4);
      cout << right << iter->second;
      aggregateContigSize += iter->first * iter->second;
      cout << " = " << aggregateContigSize << endl;
    }
  }
}


void KmerCoverageTable::printAndUpdateContigSizeHistogram(ContigNode * cNode, int tid, bool isNodeRC)
{
  if (cNode->isStandalone())
  {

    if (minValidContigLength > 0 && cNode->getLength() >= minValidContigLength)
    {
      pthread_mutex_lock(&read_file_m);
      ctgout << *cNode << endl;
      if (display_fasta == true && isNodeRC == false)
      {
        cNode->fastaDisplay(fastaout);
      }
      if (createContigSizeHistogram == true && isNodeRC == false)
      {
        uint32_t length = cNode->getLength();;
        if (contigSizeHistogram.find(length) == contigSizeHistogram.end())
        {
          contigSizeHistogram[length] = 1;
        }
        else
        {
          contigSizeHistogram[length]++;
        }
      }
      pthread_mutex_unlock(&read_file_m);
    }
    if (tight_memory_management == true)
    {
      delete cNode;
    }
    else
    {
      standaloneContigs[tid].insert(cNode);
    }
  }
  else
  {
    activeContigs[tid].insert(cNode);
  }
}


void * linkKmersParallel(void * void_args)
{
  KmerHash::iterator it;

  kmer_thread_args * args = (kmer_thread_args *) void_args;
  KmerCoverageTable * table = args->table;
  uint64_t         cLinks = 0;
  const int dotScaling    = table->dotScaling;
  int32_t localHashmapIdx = 0;
  Kmer left, right, leftRC, rightRC;

  // link kmers
  while (true)
  {
    pthread_mutex_lock(&hashmaps_m);
    localHashmapIdx = table->countHashmaps;
    table->countHashmaps++;
    pthread_mutex_unlock(&hashmaps_m);
    if (localHashmapIdx >= table->numHashmaps)
    {
      break;
    }
    KmerHash & hashmap = table->hashmaps[localHashmapIdx];
    list<KmerNode *> & graph = table->graphs[localHashmapIdx];

    for (it = hashmap.begin(); it != hashmap.end(); ++it)
    {
      KmerNode * kmerNode   = it->second.nodes;
      KmerNode * kmerRCNode = it->second.nodes + 1;
      assert(kmerNode != NULL && kmerRCNode);
      Kmer *     kmer       = it->first;
      //Kmer *     kmerRC     = kmerRCNode->getKmer();
      bool foundHasBranch   = kmerNode->hasBranch();

      for (int j = 0; j < 4; j++)
      {
        right = *kmer;
        right.pushFromLSB((nucleotide_type)j);
        rightRC = *kmer;
        rightRC.reverseComplement();
        rightRC.pushFromMSB(Kmer::ComplementaryNucleotide((nucleotide_type)j));

        bool isRightRCSmall = (rightRC < right);
        KmerNode * rightNC = (isRightRCSmall) ? table->find(&rightRC) : table->find(&right);
        if (rightNC)
        {
          KmerNode * rightNode   = rightNC;
          KmerNode * rightRCNode = rightNC + 1;
          cLinks++;
          if (right == *(rightNode->getKmer()))
          {
            kmerNode->setNext(rightNode, (nucleotide_type)j);
            kmerRCNode->setPrev(rightRCNode, Kmer::ComplementaryNucleotide((nucleotide_type)j));
          }
          else
          {
            kmerNode->setNext(rightRCNode, (nucleotide_type)j);
            kmerRCNode->setPrev(rightNode, Kmer::ComplementaryNucleotide((nucleotide_type)j));
          }
          if (foundHasBranch == false && kmerNode->hasBranch() == true)
          {
            foundHasBranch = true;
            graph.push_back(kmerNode);
            graph.push_back(kmerRCNode);
          }
        }

        left  = *kmer;
        left.pushFromMSB((nucleotide_type)j);
        leftRC = *kmer;
        leftRC.reverseComplement();
        leftRC.pushFromLSB(Kmer::ComplementaryNucleotide((nucleotide_type)j));

        bool isLeftRCSmall = (leftRC < left);
        KmerNode * leftNC = (isLeftRCSmall) ? table->find(&leftRC) : table->find(&left);
        if (leftNC)
        {
          KmerNode * leftNode   = leftNC;
          KmerNode * leftRCNode = leftNC + 1;
          cLinks++;
          if (left == *(leftNode->getKmer()))
          {
            kmerNode->setPrev(leftNode, (nucleotide_type)j);
            kmerRCNode->setNext(leftRCNode, Kmer::ComplementaryNucleotide((nucleotide_type)j));
          }
          else
          {
            kmerNode->setPrev(leftRCNode, (nucleotide_type)j);
            kmerRCNode->setNext(leftNode, Kmer::ComplementaryNucleotide((nucleotide_type)j));
          }
          if (foundHasBranch == false && kmerNode->hasBranch() == true)
          {
            foundHasBranch = true;
            graph.push_back(kmerNode);
            graph.push_back(kmerRCNode);
          }
        }
      }
    }
    if (localHashmapIdx % dotScaling == dotScaling - 1)
    {
      pthread_mutex_lock(&read_file_m);
      cout << "." << flush;
      table->nDots++;
      if (table->nDots % table->dotPerLine == 0)
      {
        cout << " processed " << table->nDots * dotScaling << " maps" << endl << flush;
      }
      pthread_mutex_unlock(&read_file_m);
    }
  }

  pthread_mutex_lock(&read_file_m);
  table->countLinks += cLinks;
  pthread_mutex_unlock(&read_file_m);
  return (void *) NULL;
}


void * populateKmerNodesParallel(void * void_args)
{
  kmer_thread_args * args = (kmer_thread_args *) void_args;
  KmerCoverageTable * table   = args->table;
  uint64_t minRepeatCount = table->minRepeatCount;
  uint64_t maxErrorCount  = table->maxErrorCount;
  const int dotScaling    = table->dotScaling;
  int32_t  tid            = args->tid;

  int32_t localHashmapIdx = 0;
  while (true)
  {
    pthread_mutex_lock(&hashmaps_m);
    localHashmapIdx = table->countHashmaps;
    table->countHashmaps++;
    pthread_mutex_unlock(&hashmaps_m);
    if (localHashmapIdx >= table->numHashmaps)
    {
      break;
    }
    KmerHash & hashmap = table->hashmaps[localHashmapIdx];
    KmerHash::iterator it = hashmap.begin();

    while (it != hashmap.end())
    {
      assert(it->first != NULL && (it->second).count > 0);
      if ((it->second).count >= minRepeatCount || (it->second).count <= maxErrorCount)
      {
        if (table->display_repeats == true && (it->second).count >= minRepeatCount)
        {
          pthread_mutex_lock(&read_file_m);
          table->ctgout << *(it->first) << ", hitcount= " << (it->second).count << endl;
          pthread_mutex_unlock(&read_file_m);
        }
        //delete it->first;
        hashmap.erase(it++);
      }
      else
      {
        /*KmerNode * newnode   = new KmerNode(it->first, 0, (it->second).count, false);
        KmerNode * newnodeRC = new KmerNode(it->first, 0, (it->second).count, true);
        it->second.nodes.kmerNode   = newnode;
        it->second.nodes.kmerNodeRC = newnodeRC;
        */uint64_t kmerNodePair = table->blockAllocators[tid]->allocate(sizeof(KmerNode)*2);
        uint64_t count = (it->second).count;
        it->second.nodes = new(kmerNodePair) KmerNode(it->first, 0, count, false);
        new(kmerNodePair+sizeof(KmerNode)) KmerNode(it->first, 0, count, true);
        ++it;
      }
    }
    if (localHashmapIdx % dotScaling == dotScaling - 1)
    {
      pthread_mutex_lock(&read_file_m);
      cout << "." << flush;
      table->nDots++;
      if (table->nDots % table->dotPerLine == 0)
      {
        cout << " processed " << table->nDots * dotScaling << " maps" << endl << flush;
      }
      pthread_mutex_unlock(&read_file_m);
    }
  }
  return (void *) NULL;
}


void * processReadsParallel(void * void_args)
{
  kmer_thread_args * args   = (kmer_thread_args *) void_args;
  KmerCoverageTable * table = args->table;
  int32_t         tid       = args->tid;
  Reader        * reader    = args->reader;
  int current               = 0;
  SequenceBuffer * buf      = table->seqbuf[tid];
  Kmer kmer, kmerRC;

  bamt.initiate();
  pthread_barrier_wait(&global_barrier);

  while (true)
  {
    pthread_mutex_lock(&read_file_m);
    current = reader->readToBuffer(buf);
    table->countReads     += current;
    table->countReadsCurr += current;
    cout << "." << flush;
    table->nDots++;
    if (table->nDots % table->dotPerLine == 0)
    {
      //cout << " read " << table->countReads << " reads, hashsize=" << table->getTableSize() << ", kmers=" << table->getKmerNumber() << endl << flush;
      cout << " read " << table->countReads << " reads" << endl << flush;
    }
    if (table->per_size_error_removal == true && 
        ((table->first_error_removal == true && table->countReadsCurr >= table->error_removal_block_size_first) ||
         (table->first_error_removal == false && table->countReadsCurr >= table->error_removal_block_size)))
    {
      cout << endl << *table << endl;
      table->countReadsCurr -= (table->first_error_removal == true ? 
                                table->error_removal_block_size_first :
                                table->error_removal_block_size);
      table->first_error_removal = false;
      for (int i = 0; i < table->numHashmaps; i++)
      {
        pthread_mutex_lock(&(table->hashmap_m[i]));
      }
      bamt.clearBlocks(false);
      bamt.initiate();
      table->hashmapstmp  = new KmerHashTmp * [table->numHashmaps];
      for (int i = 0; i < table->numHashmaps; i++)
      {
        uint64_t pointer = bamt.allocate(sizeof(KmerHashTmp));
        table->hashmapstmp[i] = new((void *)pointer) KmerHashTmp();
        pthread_mutex_unlock(&(table->hashmap_m[i]));
      }
    }
    pthread_mutex_unlock(&read_file_m);

    if (current <= 0)
    {
      return (void *) NULL;
    }
    int lineCount = buf->getLineCount();

    for (int i = 0; i < lineCount; i++)
    {
      const char * sequence = buf->getLinePointer(i);
      uint32_t len = strlen(sequence);
      int      lastInvalid = table->kmersize;

      if (len < table->kmersize)
      {
        continue;  // no kmer can be generated
      }

      // Initialize kmer
      for (uint32_t j = 0; j < len; j++)
      {
        lastInvalid = (Kmer::isValidNucleotide(sequence[j]) == false) ? table->kmersize : lastInvalid - 1;
        kmer.pushFromLSB(Kmer::charToNucleotide(sequence[j]));
        kmerRC.pushFromMSB(Kmer::charToComplementaryNucleotide(sequence[j]));

        if (lastInvalid <= 0)
        {
          if (table->per_file_error_removal == true || table->per_size_error_removal == true)
          {
            table->addParallelPrune(&kmer, &kmerRC, tid);
          }
          else
          {
            table->addParallel(&kmer, &kmerRC, tid);
          }
        }
      }
    }
    table->seqbuf[tid]->clear();
  }

  return (void *) NULL;
}


void * initializeContigsParallel(void * void_args)
{
  kmer_thread_args * args = (kmer_thread_args *) void_args;
  int32_t         tid     = args->tid;
  KmerCoverageTable * table   = args->table;
  const int dotScaling    = table->dotScaling;

  int32_t localHashmapIdx = 0;
  while (true)
  {
    pthread_mutex_lock(&hashmaps_m);
    localHashmapIdx = table->countHashmaps;
    table->countHashmaps++;
    pthread_mutex_unlock(&hashmaps_m);
    if (localHashmapIdx >= table->numHashmaps)
    {
      break;
    }
    KmerHash & hashmap = table->hashmaps[localHashmapIdx];
    KmerHash::iterator iter = hashmap.begin();

    while (iter != hashmap.end())
    {
      KmerNode * pnode = (iter->second).nodes;

      if (pnode->hasBranch() == false)
      {
        // a node with a branch is a graph node
        if (pnode->valid == false)
        {
          // kmer is moved to the ContigNode already. simply delete it
        }
        else if (pnode->prev.node == NULL || pnode->prev.node->hasBranch() == true)
        {
          uint64_t pointer = table->blockAllocators[tid]->allocate(sizeof(ContigNode));
          ContigNode * contigNode = new(pointer) ContigNode(pnode, 0);
          table->printAndUpdateContigSizeHistogram(contigNode, tid);
        }
      }

      pnode = (iter->second).nodes + 1;

      if (pnode->hasBranch() == true)
      {
        ++iter;
      }
      else
      {
        if (pnode->valid == false)
        {
          // kmer is moved to the ContigNode already. simply delete it
        }
        else if (pnode->prev.node == NULL || pnode->prev.node->hasBranch() == true)
        {
          uint64_t pointer = table->blockAllocators[tid]->allocate(sizeof(ContigNode));
          ContigNode * contigNode = new(pointer) ContigNode(pnode, 0);
          table->printAndUpdateContigSizeHistogram(contigNode, tid, true);
        }
        ++iter;
      }
    }
    if (localHashmapIdx % dotScaling == dotScaling - 1)
    {
      pthread_mutex_lock(&read_file_m);
      cout << "." << flush;
      table->nDots++;
      if (table->nDots % table->dotPerLine == 0)
      {
        cout << " processed " << table->nDots * dotScaling << " maps" << endl << flush;
      }
      pthread_mutex_unlock(&read_file_m);
    }
  }

  return (void *) NULL;
}


void * initializeContigsParallelTight(void * void_args)
{
  kmer_thread_args * args = (kmer_thread_args *) void_args;
  int32_t         tid     = args->tid;
  KmerCoverageTable * table   = args->table;
  const int dotScaling    = table->dotScaling;

  int32_t localHashmapIdx = 0;
  while (true)
  {
    pthread_mutex_lock(&hashmaps_m);
    localHashmapIdx = table->countHashmaps;
    table->countHashmaps++;
    pthread_mutex_unlock(&hashmaps_m);
    if (localHashmapIdx >= table->numHashmaps)
    {
      break;
    }
    KmerHash & hashmap = table->hashmaps[localHashmapIdx];
    KmerHash::iterator iter = hashmap.begin();

    while (iter != hashmap.end())
    {
      KmerNode * pnode = (iter->second).nodes;

      if (pnode->hasBranch() == false)
      {
        // a node with a branch is a graph node
        if (pnode->valid == false)
        {
          // kmer is moved to the ContigNode already. simply delete it
          delete pnode;       
        }
        else if (pnode->prev.node == NULL || pnode->prev.node->hasBranch() == true)
        {
          uint64_t pointer = table->blockAllocators[tid]->allocate(sizeof(ContigNode));
          ContigNode * contigNode = new(pointer) ContigNode(pnode, 0);
          table->printAndUpdateContigSizeHistogram(contigNode, tid);
          delete pnode;       
        }
      }

      pnode = (iter->second).nodes + 1;

      if (pnode->hasBranch() == true)
      {
        ++iter;
      }
      else
      {
        if (pnode->valid == false)
        {
          // kmer is moved to the ContigNode already. simply delete it
          delete pnode;       
        }
        else if (pnode->prev.node == NULL || pnode->prev.node->hasBranch() == true)
        {
          uint64_t pointer = table->blockAllocators[tid]->allocate(sizeof(ContigNode));
          ContigNode * contigNode = new(pointer) ContigNode(pnode, 0);
          table->printAndUpdateContigSizeHistogram(contigNode, tid, true);
          delete pnode;       
        }
        ++iter;
      }
    }
    hashmap.clear();
    if (localHashmapIdx % dotScaling == dotScaling - 1)
    {
      pthread_mutex_lock(&read_file_m);
      cout << "." << flush;
      table->nDots++;
      if (table->nDots % table->dotPerLine == 0)
      {
        cout << " processed " << table->nDots * dotScaling << " maps" << endl << flush;
      }
      pthread_mutex_unlock(&read_file_m);
    }
  }

  return (void *) NULL;
}


void KmerCoverageTable::runInteractiveMode()
{
  if (interactive_mode == false)
  {
    return;
  }

  while (true)
  {
    cout << "type one of followings:" << endl;
    cout << " - kmer node number starting with kn"   << endl;
    cout << " - contig node number starting with cn" << endl;
    cout << " - kmer value (capital letters)"        << endl;
    cout << " - q or quit for exiting the program"   << endl;

    string input;
    getline(cin, input);

    if (input == "q" || input == "quit")
    {
      exit(1);
    }
    if (input.size() < 2)
    {
      continue;
    }
    else if (input.substr(0, 2) == "cn")
    {
      istringstream cnNumStr(input.substr(2));
      uint64_t      cnNum = 0;
      bool          found = false;

      if (cnNumStr.good() == true)
      {
        cnNumStr >> cnNum;
      }
      else
      {
        continue;
      }
      for (int i = 0; i < numThreads; i++)
      {
        ContigHash & chash = activeContigs[i];
        for (ContigHash::iterator iter = chash.begin(); iter != chash.end(); ++iter)
        {
          if (cnNum == (*iter)->id)
          {
            cout << *(*iter) << endl;
            found = true;
            i     = numThreads;
            break;
          }
        }
      }
      if (found == true)
      {
        continue;
      }
      for (int i = 0; i < numThreads; i++)
      {
        ContigHash & chash = standaloneContigs[i];
        for (ContigHash::iterator iter = chash.begin(); iter != chash.end(); ++iter)
        {
          if (cnNum == (*iter)->id)
          {
            cout << *(*iter) << endl;
            found = true;
            i     = numThreads;
            break;
          }
        }
      }
    }
    else if (input.substr(0, 2) == "kn")
    {
      cout << "kn is not yet supported" << endl;
    }
    else if (input.size() == kmersize)
    {
      Kmer kmer;
      Kmer kmerRC;

      for (uint32_t i = 0; i < kmersize; i++)
      {
        kmer.pushFromLSB(Kmer::charToNucleotide(input[i]));
        kmerRC.pushFromMSB(Kmer::charToComplementaryNucleotide(input[i]));
      }

      if (kmerRC < kmer)
      {
        kmer = kmerRC;
      }

      int hashValue         = kmer.getLSW() % numHashmaps;
      KmerHash & hashmap    = hashmaps[hashValue];
      KmerHash::iterator it = hashmap.find(&kmer);

      if (it != hashmap.end())
      {
        cout << *((it->second).nodes)   << endl;
        cout << *((it->second).nodes + 1) << endl;
      }
    }
  }
}

